/*
 heat-spectral lib. solves time-dependent PDEs using SEM-NI method.
 Copyright (C) 2009, Altieri Michele and Stefanoni Andrea
 Copyright (C) 2008, Proverbio Alessandro and Radice David

 This program is free software: you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
 the Free Software Foundation, either version 3 of the License, or
 (at your option) any later version.
 
 This program is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 GNU General Public License for more details.
 
 You should have received a copy of the GNU General Public License
 along with this program. If not, see <http://www.gnu.org/licenses/>.
*/

/**
   \file mesh.hpp
   \brief This file contains the Mesh class

   The Mesh class contains data and methods for the rappresentation
   and the use of a two dimensional SEM mash over a rectangular domain
*/

# ifndef MESH_HPP
# define MESH_HPP 1

# include<cmath>
# include<cassert>
# include<fstream>
# include<valarray>
# include<vector>

# include "point.hpp"
# include "integrator.hpp"
# include "map2d.hpp"
# include "functions.hpp"
# include "solution.hpp"
# include "typedefs.hpp"

namespace spectral {

/**
    \class Mesh
    \brief Mesh Class

    This class contains data and methods for the rappresentation and
    the use of a two dimensional SEM mash over a rectangular domain,
*/

/**
    \fn Mesh::Mesh(int nx,int ny,int nnodes, int nmax,Point<double> const & A,
                    Point<double> const & B,double legDelta,double legEps, int order=1 )
    \brief Mesh Constructor
    \param nx Number of elements on the x-axis
    \param ny Number of elements on the y-axis
    \param nnodes Number of nodes for the quadrature in every element
    \param nmax Maximum number of iteration for the Newton-Horner method
    \param Point<double> A Bottom-Left point of the dominium
    \param Point<double> B Top-Right point of the dominium
    \param legDelta Minimum increment in the newton formula used to compute the GLL nodes
    \param legEps Critical resiudual in the newton formula used to compute the GLL nodes
    \param order Is the order of the problem in the week form. For example it's 1 for the standard advection-diffusion problem.

    This is the constructor.
*/

/**
    \fn Mesh::Mesh(Mesh const & Me)
    \brief Mesh copy constructor

    It justs make a copy of a given Mesh
*/

/**
    \fn Mesh::Mesh(Mesh const & Me,Point<int> const & A, Point<int> const & B)
    \brief Sub-Mesh constructor
    \param A The coordinates of the bottom-left point of
           the new mesh. Coordinates are the integer indexes on the lattice of
           the mesh
    \param B The coordinates of the top-right point of
           the new mesh. Coordinates are the integer indexes on the lattice of
           the mesh


    This will build a sub-mesh of a given mesh.
*/

/**
    \fn Mesh::Mesh(Mesh const & Me,Point<double> const & A, Point<double> const & B)
    \brief Sub-Mesh constructor
    \param A The coordinates of the bottom-left point of
           the new mesh.
    \param B The coordinates of the top-right point of
           the new mesh.

    This will build a sub-mesh of a given mash. The real coordinates of the new
    subgrid are taken as the best approximation of the given coordinates on the
    mesh lattice.
*/

/**
    \fn Mesh::~Mesh()
    \brief Destructor

*/

/**
    \fn double Mesh::prodL2(T const & t, F const & f) const
    \brief Computes the discrete scalar product

    This will calculate the discrete scalar product as:
    \f[ \sum_{j=1}^M \big( t , f \big)_{\Delta,j} \f]
    where \e j runs on the number of elements and \f$(t,f)_{\Delta,j}\f$ is the
    GLL approximation of the integral of the product over the j-sim element.
*/

/**
    \fn double Mesh::prodL2(T const & t, F const & f, int elemnum) const
    \brief Computes the discrete scalar product on the single element

    This will calculate the discrete scalar product of t and f
    on the \e elemnum element of the SEM mesh. t and f must be two functors
    with the operator() method overloaded with signature:
    \code double operator()(double,double) \endcode
*/

/**
    \fn double Mesh::prodL2(T1 const & t1, T2 const & t2, T3 const & t3, int elemnum) const
    \brief Compute the discrete scalar product on the single element

    This will compute the discrete scalar product of t1,t2 and t3
    on the \e elemnum element of the SEM mesh. t1, t2 and t3 must be
    functors with the operator() method overloaded with signature:
    \code double operator()(double,double) \endcode
*/

/**
    \fn double Mesh::prodL2border(T const & t, F const & f, int elemnum, int border) const
    \brief Computes the discrete scalar product on a single border of a specific element

    This will calculate the discrete scalar product of t and f on the
    \e border border of the \e elemnum element. t and f must be two functors with
    the operator() method overloaded with signature:
    \code double operator()(double,double) \endcode
*/

/**
    \fn double Mesh::integrate(std::valarray<double> const & v,
                               int elemnum) const
    \brief Integrates \e v over the element \e elemnum

    This is to be used when it's necessary to integrate a function
    from its interpolation on the grid.

    The functor approach is to be prefered everywhere else!
    (This method is much slower and less flexible).
*/

/**
    \fn double Mesh::numNorm(std::valarray<double> const & v) const
    \brief Numeric norm of v

    The numeric norm is computed as:
    \f[
    \| v \|^2_{\Delta} = \sum_{e=0}^{N-1} \sum_{i=0}^{K-1} \alpha_{e,i} v_{e,i}
    \f]
    if \e N is the number of elements in the mesh and \e K is the
    number of nodes in the number of GLL nodes in a single element.
*/

/**
    \fn int Mesh::getnx() const
    \brief Returns the number of elements along the x-axis
*/

/**
    \fn int Mesh::getny() const
    \brief Returns the number of elements along the y-axis

*/

/**
    \fn int Mesh::getnelem() const
    \brief Returns the total number of elements in the mesh

*/

/**
    \fn int Mesh::getnnodes() const
    \brief Returns the number of GLL nodes

    This will return the number of nodes for the GLL quadrature
    along each axis. So the total number of nodes in each element
    is \e (Mesh::getnnodes())*(Mesh::getnnodes()).
*/

/**
    \fn int Mesh::getProblemDim() const
    \brief Returns the number of degrees of freedom on the whole mesh

*/

/**
    \fn int Mesh::getnbc() const
    \brief Returns the number of nodes on the boundary of the domain
*/

/**
    \fn int Mesh::getnNodesPerEl() const
    \brief Returns the number of nodes in a single spectral element

*/

/**
    \fn int Mesh::getnNodesPerBord(int border) const
    \brief Returns the number of nodes on the \e border boundary.

   The convenction for the border is the same of the
   Mesh::getBordElements function.
*/

/**
    \fn void Mesh::getElement(int elemnum,Point<double> & BL,Point<double> & TR) const
    \brief Returns the bottom-left and the top-right points of the given element

*/

/**
 * \fn int Mesh::getGlobalNum(int elemnum,int idxnode ) const
 * \brief Returns the global numeration of the given node of the given element
 */

/**
 * \fn int Mesh::getGlobalNum(int elemnum,int idx,int idy) const
 * \brief Returns the global numeration of the given node of the given element
 */

/**
    \fn int Mesh::getNodeNumber(int idx,int idy) const
    \brief Returns the global node numeration

    If \e idx and \e idy are the cartesian coordinates on the nodes
    lattice, this method will return the global numeration of the
    given node
*/

/**
    \fn void Mesh::getGlobalNum(int elemnum,
                                std::vector<int> & nodes) const
    \brief Gets the global numeration for the given element

    This method will store the global index values for the nodes of the
    given element into \e nodes. The chosen order is from left to right
    and from the bottom to the top so that if \e i is a GLL node in the
    local numeration \e nodes[i] will contain the number of the node in
    the global numeration.
*/

/**
   \fn void Mesh::getBordElements(int border, std::vector<int> & elements) const
   \brief Elements with at least one vertex on the border

   This method will store the list of the elements with at least one node
   on the border \e border. The convenction for the \e border is:
   \code
                  2
      ------------------------
      |                      |
      |                      |
    3 |                      |  1
      |                      |
      |______________________|

                  0
   \endcode
 */

/**
   \fn void Mesh::getBordNodes(int border, std::vector<int> & nodes) const
   \brief Vertexes on the border of an element

   This method will store the list of nodes on the border \e border. The
   convenction for the border is the same of the Mesh::getBordElements function.
   The vertex number are in the \b local numeration.
 */

/**
    \fn void Mesh::whichBord(int node,bool isonborder[4]) const
    \brief Which edges owns the given point

    \e isonborder[i] will be \e true if the node \e node is on border \e i
*/

/**
    \fn bool Mesh::isPhysicalVertex(int node) const
    \brief True if the given node is the vertex of the physical domain

*/

/**
    \fn int Mesh::getBordNum(int border, int i) const
    \brief Returns the global numeration of the node i on the border \e border

    \e i must be a node on the border, numbered with a counter-clock-wise
    convenction. \e i must go from 0 to the number of nodes along the border.
    \see{Mesh::getnNodesPerBord}
*/

/**
    \fn void Mesh::getPhysicalVertex(int node,int & bord1,int & bord2) const
    \brief Gives the edges containing the given vertex

    \e node \b must \b be a physical vertex otherwise this method could
    go in an endless loop.
*/

/**
   \fn void Mesh::getGlobBordNodes(int border, int elemnum, std::vector<int> & nodes) const
   \brief Vertexes on the border of an element

   This method will store the list of nodes on the border \e border. The
   convenction for the border is the same of the Mesh::getBordElements function.
   The vertex number are in the \b global numeration.

   Please mind that this method will assume that the given element is on the
   given border, no check is made to control this.
*/

/**
    \fn Point<double> Mesh::getNodePosition(int node) const
    \brief Coordinates of the given node in the physical domain

    \e node must be a valid node in the \b global numeration.
*/

/**
    \fn Point<double> Mesh::getNodePosition(int pdx,int pdy) const
    \brief Coordinates of the given node in the physical domain

    \e pdx and \e pdy must be the coordinates of the wanted node
    in the nodes lattice.
*/

/**
    \fn Point<double> Mesh::getRefNodePosition(int node) const
    \brief Coordinates of the given node in the reference element

    \e node must be a valid node in the \b global numeration.
*/

/**
    \fn Point<double> Mesh::getRefNodePosition(int pdx,int pdy) const
    \brief Coordinates of the given node in the reference element

    \e pdx and \e pdy must be the coordinates of the wanted node
    in the nodes lattice.
*/

/**
    \fn Point<int> Mesh::getVertexPosition(int edx, int edy) const
    \brief Coordinates of the given vertex in the nodes lattice

    This is used to convert the coordinates of a given point in the
    large grid of the spectral element to the coordinates in the
    fine grid of the GLL nodes.
*/

/**
    \fn int whereami(Point<double> const & P) const
    \brief Element number which contains the given point

    The point \e P must be a valid point in the physical domain.
*/

/**
    \fn Point<int> Mesh::gridPosition(int node) const
    \brief Coordinates of the given node in the nodes lattice
*/

/**
    \fn Point<int> Mesh::gridPosition(Point<double> const & P) const
    \brief Coordinates of a point on the grid

    This will return the coordinates of a point on the lattice grid of the SEM mesh.
    The indexes of the best approximating vertices is returned. Please notice
    that this will NOT return the best approximate node, but the best aproximating
    vertex.
*/


/**
    \fn int Mesh::whereami(Point<double> const & P) const
    \brief Which element contains the given point

    This will return the number of the element containing the given point.
*/

/**
    \fn inline Functor<T,Map2d<direct> > Mesh::restriction(T const & f,
                                                           int elemnum) const
    \brief Returns the restriction of a given function on a given element

    This is to be used to build a wrapper to a given function before the
    calculation of the integral. The returned object will be a composition
    of the function and the affine transformation. \see{Functor}.
*/

/**
    \fn inline Solution Mesh::restriction(std::valarray<double> const & u,
                                          int elemnum) const
    \brief Returns the restriction of the given solution vector on a given element

    This is used to interpolate the solution on a given element in order
    to further evaluate it or compute numerical integration or numeric
    differentiation.
*/

/**
    \fn Functor<T, Map1d<direct> > Mesh::trace(T const & f,
                                               int border) const
    \brief Trace operator of f on the border \e border.

    This will return a Functor which is used to
    evaluate f over the border to perform numerical
    integration. The convenction for the border is:
    \code
                    2
        ------------------------
        |                      |
        |                      |
      3 |                      |  1
        |                      |
        |______________________|

                    0
    \endcode
*/

/**
    \fn Point<int> Mesh::meshDimension( void )
    \brief Returns the dimension in number of element of the mesh

    This function is useful for extract the data of the mesh. It has in output a Point<int>.
*/

/**
    \fn int Mesh::meshDegree( void )
    \brief Returns the degree of the polynomial used in the mesh

    Like meshDimension is useful for take out from the mesh their parameters.
*/

/**
	\fn TensorProd Mesh::lagrange(int i) const
	\brief Returns \f$ \hat{l}_i \f$

	This returns an element of the lagrangian basis on the reference interval
	\f$ [-1,1] \f$
*/

/**
    \fn TensorProd Mesh::lagrange(int i, int j) const
    \brief Returns \f$ \hat{l}_{ij} \f$

    This returns an element of the lagrangian basis on the reference square
    \f$ [-1,1] \times [-1,1] \f$
*/

/**
    \fn Functor<TensorProd, Map2d<inverse> > Mesh::lagrange(int i, int j,
                                                            int elemnum) const
    \brief Returns \f$ \hat{l}_{ij} \f$

    This returns the composition of the element of the lagrangian basis with
    the affine transformation from the reference square \f$ [-1,1] \times [-1,1] \f$
    to the given spectral element.
*/

/**
    \fn void Mesh::elementNumber(int n, Point<double>  & P1, Point<double> & P2) const
    \brief Computes the edges of a given spectral element

    After the execution of the function P1 will contain the bottom-left point of the
    n-sim element, while P2 will contain the top-right point.
*/

/**
    \fn inline TensorProd Mesh::Dx(TensorProd const & tp, int elemnum) const
    \brief Partial derivatives of tp with respect to x

    This should be used everytime that we want to compute the L2
    scalar product between a function and a derivative of an element
    of the basis. The derivative is computed on the fisical domain.
*/

/**
    \fn inline Solution Mesh::Dx(Solution const & sol,int elemnum) const
    \brief Partial derivatives of sol with respect to x

    This should be used everytime that we want to compute the L2
    scalar product between a function and a derivative of the interpolated
    solution. The derivative is computed on the fisical domain.
*/

/**
    \fn inline TensorProd Mesh::Dy(TensorProd const & tp, int elemnum) const
    \brief Partial derivatives of tp with respect to y

    This should be used everytime that we want to compute the L2
    scalar product between a function and a derivative of an element
    of the basis. The derivative is computed on the fisical domain.
*/

/**
    \fn inline Solution Mesh::Dy(Solution const & sol,int elemnum) const
    \brief Partial derivatives of sol with respect to y

    This should be used everytime that we want to compute the L2
    scalar product between a function and a derivative of the interpolated
    solution. The derivative is computed on the fisical domain.
*/

/**
    \fn void Mesh::save_to_file(std::valarray<double> const & vec,
                                char const * filename) const
    \brief Writes the given array in the given file

    The output file format is an ASCII formatted to be used with the
    GNUPLOT plotting function "splot". Consult the GNUPLOT manual for
    more informations about this format.
*/

/**
    \fn void Mesh::save_to_file(double const * vec,
                                char const * filename) const
    \brief Writes the given array in the given file

    The output file format is an ASCII formatted to be used with the
    GNUPLOT plotting function "splot". Consult the GNUPLOT manual for
    more informations about this format.
*/

/**
    \fn void Mesh::interpolate(T const & f,
                               std::valarray<double> & v) const
    \brief Interpolates \e f on the mesh

    \e v will contain the values of \e f on the nodes of the grid.

    This function is \b not to be used to build an array to be used
    to perform numerical integration: the functor apprach is to used
    for this.
*/

class Mesh {
    public:
        // Constructors and desctructor
        // ----------------------------
        Mesh(int nx,
             int ny,
             int nnodes,
             int nmax,
             Point<double> const & A,
             Point<double> const & B,
             double legDelta,
             double legEps,
             int order=1);

        Mesh(Mesh const & Me);

        Mesh(Mesh const & Me,
             Point<int> const & A,
             Point<int> const & B);

        Mesh(Mesh const & Me,
             Point<double> const & A,
             Point<double> const & B);

        ~Mesh();

        // L2 scalar product API
        // ---------------------
        template<class T, class F>
        double prodL2(T const & t,
                      F const & f) const;

        template<class T, class F>
        double prodL2(T const & t,
                      F const & f,
                      int elemnum) const;

        template<class T1, class T2, class T3>
        double prodL2(T1 const & t1,
                      T2 const & t2,
                      T3 const & t3,
                      int elemnum) const;

	template<class T, class F>
	double prodL2border(T const & t,
			    F const & f,
			    int elemnum,
			    int border) const;

        double integrate(std::valarray<double> const & v,
                         int elemnum) const;

        double numNorm(std::valarray<double> const & v) const;

        // Geometric methods
        // -----------------
        int getnx() const;
        int getny() const;
        int getnelem() const;
        int getnnodes( void ) const;
        int getProblemDim() const;
        int getnbc() const;
        int getnNodesPerEl() const;

        int getnNodesPerBord(int border) const;

        void getElement(int elemnum,
                        Point<double> & BL,
                        Point<double> & TR) const;

        int getGlobalNum(int elemnum,
                         int idxnode ) const;

        int getGlobalNum(int elemnum,
                         int idx,
                         int idy) const;

        void getGlobalNum(int elemnum,
                          std::vector<int> & nodes) const;

        int getNodeNumber(int idx, int idy ) const;

        void getBordElements(int border,
                             std::vector<int> & elements) const;

        void getBordNodes(int border,
                          std::vector<int> & nodes) const;

	void getDomainBordNodes(int border,
				std::vector<int> & nodes) const;

        void whichBord(int node,
                       bool isonborder[4]) const;
        bool isPhysicalVertex(int node) const;
        void getPhysicalVertex(int node,int & bord1,int & bord2) const;

        int getBordNum(int border,int i) const;

        void getGlobBordNodes(int border,
                              int elemnum,
                              std::vector<int> & nodes) const;

        Point<double> getNodePosition(int node) const;
        Point<double> getNodePosition(int pdx,int pdy) const;
        Point<double> getRefNodePosition(int node) const;
        Point<double> getRefNodePosition(int pdx,int pdy) const;

        Point<int> getVertexPosition(int edx,int edy) const;

        Point<int> gridPosition(int node) const;
        Point<int> gridPosition(Point<double> const & P) const;

        int whereami(Point<double> const & P) const;

        void elementNumber(int n,
                           Point<double> & P1,
                           Point<double> & P2) const;

        // Functors operations
        // -------------------
        template<class T>
        inline Functor<T, Map2d<direct> > restriction(T const & f,
                                                      int elemnum) const;

        inline Solution restriction(std::valarray<double> const & u,
                                    int elemnum) const;

	inline TensorProd lagrange(int i) const;

	inline TensorProd lagrange(int i, int j) const;

        inline TensorProd Dx(TensorProd const & tp,
                             int elemnum) const;

        inline Solution Dx(Solution const & sol,
                           int elemnum) const;

        inline TensorProd Dy(TensorProd const & tp,
                             int elemnum) const;

        inline Solution Dy(Solution const & sol,
                           int elemnum) const;

        inline
        Functor<TensorProd,Map2d<inverse> > lagrange(int i,
                                                     int j,
                                                     int elemnum) const;

        // Misc Functions
        // --------------
        template<class T>
        void interpolate(T const & f,
                         std::valarray<double> & v) const;

        void save_to_file(std::valarray<double> const & vec,
                          char const * filename) const;
        void save_to_file(double const * vec,
                          char const * filename) const;

        // Friends
        // -------
        friend class Solution;

    private:
        int _nx;              /**< \brief Number of elements (not vertex!) along the x-axis */
        int _ny;              /**< \brief Number of elements (not vertex!) along the y-axis */
        int _nnodes;          /**< \brief Number of nodes for the GLL quadrature */
        Point<double> _A;     /**< \brief Bottom-left point */
        Point<double> _B;     /**< \brief Top-right point */
        double _legDelta;     /**< \brief Legendre polynomial computational delta */
        double _legEps;       /**< \brief Legendre polynomial computational epsilon */

        double * _vertexX;    /**< \brief Vertexes on the x-axis */
        double * _vertexY;    /**< \brief Vertexes on the y-axis */
        Integrator _integr;   /**< \brief Integrator class */
};

// ====================================================
// End of proper header
// ====================================================

template<class T, class F>
double Mesh::prodL2(T const & t, F const & f) const {
    Product<T,F> prod(t,f);
    double sum=0;
    Point<double> C(double(0),double(0));
    Point<double> D(double(0),double(0));
    for(int i=0;i<_nx;++i) {
        for(int j=0;j<_ny;++j) {
            C[0]=_vertexX[i];
            C[1]=_vertexY[i];
            D[0]=_vertexX[i+1];
            D[1]=_vertexY[i+1];

            sum+=_integr(prod,C,D);
        }
    }
    return sum;
}

template<class T, class F>
double Mesh::prodL2border(T const & t,
			  F const & f,
			  int elemnum,
			  int border) const {
    Product<T,F> prod(t,f);
    Point<double> C;
    Point<double> D;
    elementNumber(elemnum,C,D);

    double xi(0);
    double a(0);
    double b(0);

    switch(border) {
	case 0: xi = C[1];
		a = C[0];
		b = D[0];
		break;
	case 1: xi = D[0];
		a = C[1];
		b = D[1];
		break;
	case 2: xi = D[1];
		a = C[0];
		b = D[0];
		break;
	case 3: xi = C[0];
		a = C[1];
		b = D[1];
		break;
	default:
		exit(1);
    }

    return _integr(prod,a,b,xi,border);
}

template<class T, class F>
double Mesh::prodL2(T const & t,
                    F const & f,
                    int elemnum) const {
    Product<T,F> prod(t,f);
    Point<double> C;
    Point<double> D;
    elementNumber(elemnum,C,D);
    return _integr(prod,C,D);
}

template<class T1, class T2, class T3>
double Mesh::prodL2(T1 const & t1,
                    T2 const & t2,
                    T3 const & t3,
                    int elemnum) const {
    Product<T2,T3> prod(t2,t3);
    return prodL2(t1,prod,elemnum);
}

template<class T>
inline Functor<T, Map2d<direct> > Mesh::restriction(T const & f,
                                                    int elemnum) const {
    Point<double> C;
    Point<double> D;
    elementNumber(elemnum,C,D);
    Map2d<direct> phi(C,D);
    return Functor<T,Map2d<direct> >(f,phi);
}

inline Solution Mesh::restriction(std::valarray<double> const & u,
                                  int elemnum) const {
    return Solution(&u[0],elemnum,*this);
}

inline TensorProd Mesh::lagrange(int i) const {
	return _integr.lagrange(i);
}

inline TensorProd Mesh::lagrange(int i, int j) const {
    return _integr.lagrange(i,j);
}

inline TensorProd Mesh::Dx(TensorProd const & tp, int elemnum) const {
    Point<double> A;
    Point<double> B;
    elementNumber(elemnum,A,B);
    double lengthX=B.get(0)-A.get(0);
    assert(lengthX>0);
    return _integr.Dx(tp,lengthX);
}

inline Solution Mesh::Dx(Solution const & sol, int elemnum) const {
    Solution dsol(sol);
    dsol.Dx(elemnum);
    return dsol;
}

inline TensorProd Mesh::Dy(TensorProd const & tp, int elemnum) const {
    Point<double> A;
    Point<double> B;
    elementNumber(elemnum,A,B);
    double lengthY=B.get(1)-A.get(1);
    assert(lengthY>0);
    return _integr.Dy(tp,lengthY);
}

inline Solution Mesh::Dy(Solution const & sol, int elemnum) const {
    Solution dsol(sol);
    dsol.Dy(elemnum);
    return dsol;
}

inline Functor<TensorProd, Map2d<inverse> > Mesh::lagrange
                                (int i, int j, int elemnum) const {
    Point<double> C;
    Point<double> D;
    elementNumber(elemnum,C,D);
    Map2d<inverse> phi(C,D);
    return Functor<TensorProd, Map2d<inverse> >(lagrange(i,j),phi);
}

template<class T>
void Mesh::interpolate(T const & f,
                       std::valarray<double> & v) const {
    Point<double> P;
    v.resize(getProblemDim());

    int _vertex;

    // Number of nodes along the x-axis
    int nnx=_nx*_nnodes-_nx+1;
    // Number of nodes along the y-axis
    int nny=_ny*_nnodes-_ny+1;

    for(int idx=0;idx<nnx;++idx) {
        for(int idy=0;idy<nny;++idy) {
             _vertex=idx+idy*nnx;
             P=getNodePosition(idx,idy);
             v[_vertex]=f(P[0],P[1]);
        }
    }
}

} // namespace

#endif
